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Abstract. The double-differential production cross-section of positive pions, (Pa^^ /dpdSl, measured in the HARP 
experiment is presented. The incident particles are 8.9 GeV/c protons directed onto a beryllium target with a thickness 
of 5% of a nuclear interaction length. The measured cross-section has a direct impact on the prediction of neutrino 
fluxes for the MiniBooNE and SciBooNE experiments at Fermilab. After cuts, 13 million protons on target produced 
about 96,000 reconstructed secondary tracks which were used in this analysis. Cross-section results are presented in the 
kinematic range 0.75 GeV/c < Ptv < 6.5 GeV/c and 30 mrad < < 210 mrad in the laboratory frame. 
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1 Introduction 

The HARP experiment was designed to make measurements 
of hadron yields from a large range of nuclear targets and for in- 
cident particle momenta from 1.5 GeV/c - 15 GeV/c. Among 
its primary goals were to contribute to the detailed understand- 
ing of neutrino beams of several experiments, including: 

- The K2K experiment, which has recently published its final 
results ||2| confirming the evidence of atmospheric oscilla- 
tions observed by Super-Kamiokande |3 1. 

- The MiniBooNE experiment |i4J, which recently excluded 
||5l two neutrino appearance-only oscillations as an expla- 
nation of the LSND anomaly (6), in the hypothesis that 
the oscillations of neutrinos and antineutrinos are the same. 
The MiniBooNE detector will also be used to measure neu- 
trino interaction cross-sections for which an absolute pre- 
diction of neutrino fluxes becomes of particular importance. 

- The SciBooNE experiment |9|, which will take data in the 
same neutrino beam used by MiniBooNE in order to per- 
form a precision measurement of neutrino cross-sections in 
the energy region around 1 GeV. 

The calculation of the flux and relative neutrino composi- 
tion of a neutrino beam requires a precise measurement of the 
interaction cross-section between the beam particles and the 
target material. In the case of the K2K and the MiniBooNE and 
SciBooNE experiments, the dominant component of the beam 
(muon neutrinos) comes from the decay of positive pions pro- 
duced in the collisions of incident protons on a nuclear target. 
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To compute the t'^ flux one needs a 47r parameterization of the 
differential cross section, cPa'^^ / dpdfl, which, in order to be 
reliable, must be based on a wide-acceptance, precise measure- 
ment. The physics program of the HARP experiment includes 
the measurement of these cross-sections. 

An earlier publication reported measurements of the 7r+ 
cross-sections from an aluminum target at 12.9 GeV/c lITOl . 
This corresponds to the energies of the KEK PS and the tar- 
get material used by the K2K experiment. The K2K oscillation 
result relies on both the measurement of an overall deficit of 
muon neutrino interactions and on the measurement of an en- 
ergy spectrum deformation observed at the Super-Kamiokande 
far detector compared to the no-oscillations expectations. In- 
troducing the HARP experimental input in the K2K oscilla- 
tion analysis has been particularly beneficial in reducing the 
systematic uncertainty in the overall number of muon neutrino 
interactions expected; the near-to-far flux extrapolation contri- 
bution to this uncertainty was reduced from 5.1% [1| to 2.9% 
L2J. 

Our next goal is to contribute to the understanding of the 
MiniBooNE and SciBooNE neutrino fluxes. They are both pro- 
duced by the Booster Neutrino Beam at Fermilab which orig- 
inates from protons accelerated to 8.9 GeV/c by the Fermi- 
lab Booster before being collided against a beryllium target. 
As was the case for the K2K beam, an important input for the 
calculation of the resulting flux is the 7r+ production cross- 
sections from a beryllium target at 8.9 GeV/c, which will be 
presented in this paper 

The HARP experimental apparatus is effectively divided 
into two tracking and particle identification sub-systems, a small- 
angle/high-momentum detection system {0: 0-0.25 rad, p: 0.5- 
8 GeV/c) and a large-angle/low-momentum system (0: 0.35- 
2.15 rad, p: 0.1-0.8 GeV/c). Fig. [T] shows a schematic of the 
HARP detector. Five modules of the NOMAD drift chambers 
[8 1 (NDCl-5) and the dipole magnet comprise the forward spec- 
trometer; a time-of-flight wall (TOFW), Cherenkov detector 
(CHE) and electromagnetic calorimeter (ECAL) make up the 
particle identification (PID) system. The large angle tracking 
and PID system is comprised of a time projection chamber 
(TPC) and resistive plate chambers (RPCs). The relevant me- 
son production for the creation of the MiniBooNE and Sci- 
BooNE neutrino fluxes is forward (0-0.30 rad) and at large 
momenta (0.5-6 GeV/c). These ranges are best covered by the 
forward tracking system and PID detectors and so the large an- 
gle system is not used in the present analysis. 

The results reported here are based on data taken in 2002 in 
the T9 beam of the CERN PS. About 2.3 million incoming pro- 
tons were selected. After cuts, 95,897 reconstructed secondary 
tracks were used in the analysis. The absolute normalization 
of the cross-section was determined using 204,295 'minimum- 
bias' trigger events. 

The analysis used in the calculation of the p-Be 7r+ produc- 
tion cross-sections being presented here follows largely from 
that used in a previous publication of p-Al 7r+ production cross- 
sections |10|. The present analysis description, therefore, will 
focus on the differences in the analysis compared to the p-Al 
publication. 

This paper is organized as follows. In Section |2] we sum- 
marize the main changes made to the analysis since the p-Al 
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Drift Chambers 




Fig. 1. Schematic layout of the HARP spectrometer. The convention for the coordinate system is shown in the lower-right corner. The three 
most downstream (unlabeled) drift chamber modules are only partly equipped with electronics and not used for tracking. 



publication. In Section |3] we describe the calculation of the 
cross-section and define its components. The following three 
sections expand on aspects of the analysis where significant 
changes have been made since the previous publication. Sec- 
tion|4]describes event and track selection and reconstruction ef- 
ficiencies. Section|5]discusses the determination of the momen- 
tum resolution and scale in the forward spectrometer. Section 
|6]summarizes the particle identification techniques. Physics re- 
sults are presented in Section [T] Section |8] discusses the rele- 
vance of these results to neutrino experiments. Finally, a sum- 
mary is presented in Section|9] 

2 Summary of analysis changes since the 
HARP p-AI publication 

The analyses of the 12.9 GeV /c p-Al data and the 8.9 GeV /c 
p-Be data are largely the same. To avoid repetition of informa- 
tion, the reader is referred to that earlier publication for many 
details not directly discussed in the present paper. The sec- 
tions concerning the experimental apparatus, the description 
of the tracking algorithm for the forward spectrometer and the 
method of calculating the track reconstruction efficiency are 
all directly valid here. The method of particle identification has 
not changed; it is only the PID detector hit selections and there- 
fore their response functions which have been significantly im- 
proved. The most important improvements introduced in this 
analysis compared with the one presented in [10] are: 

- An improvement in the minimization performed as part 
of the tracking algorithm has eliminated the anomalous dip 
in tracking efficiency above 4 GeV/c shown in |10|. The 
tracking efficiency is now > 97% everywhere above 2 GeV/c. 
(See Sec.l43Tl. 

- Studies of HARP data other than that described here have 
enabled a validation of our Monte Carlo simulation of low- 
energy hadronic interactions in carbon. Specifically, we have 



compared low energy pH-C and tt+C cross-sections to dis- 
tributions from the Binary cascade |13 | and Bertini intra- 
nuclear cascade 1 14] hadronic interaction models used to 
simulate the secondary interactions of p, n and tt*. The 
material in the HARP forward spectrometer where tertiary 
tracks might be produced is predominantly carbon. Conse- 
quently, the systematic error on the subtraction of tertiary 
tracks has been reduced from 100% in fTU\ to 50%. (See 
Sec.gaii. 

- Analysis techniques were developed for comparing the mo- 
mentum reconstructions in data and Monte Carlo allowing 
data to be used to fine-tune the drift chamber simulation pa- 
rameters. These efforts have reduced the momentum scale 
uncertainty from 5% in ifTOl to 2% in the present analy- 
sis and provided a better understanding of the momentum 
smearing caused by the HARP spectrometer, including our 
knowledge of the non-Gaussian contributions to the resolu- 
tion function. (See Sec.|5]l. 

- New selection cuts for PID hits in TOFW and in CHE have 
resulted in much reduced backgrounds and negligible effi- 
ciency losses. Consequently, the uncertainty on the cross- 
section arising from particle identification was reduced by 
a factor of seven to 0.5% making PID now a negligible con- 
tribution to the systematic error in pion yield measurements 
at forward angles. (See Sec.|6l). 

- Improved knowledge of the proton beam targeting efficiency 
and of fully correlated contributions to track reconstruc- 
tion and particle identification efficiencies have reduced the 
overall normalization uncertainty on the pion cross-section 
measurement from 4% to 2%. 

- Significant increases in Monte Carlo production have re- 
duced uncertainties from Monte Carlo statistics and allowed 
studies to reduce certain systematics to be made. 

The statistical precision of the data, however, is noticeably 
roTse. The 8.9 GeV/c beryllium and empty-target data sets 
re both smaller than the corresponding 12.9 GeV/c sets, with 
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73% and 42% of the protons-on-target for the target and empty- 
target configurations, respectively. The statistics of the target 
sample is further reduced by the p-Be total interaction cross- 
section being roughly 50% of the p-Al total cross-section. 

In the present paper the p-Be cross-sections are presented in 
13 momentum bins from 0.75 GeV/c- 6.5 GeV/c whereas the 
p-Al cross-sections were presented in 8 bins. This new binning 
was selected to attain roughly equal statistical and systematic 
uncertainties - on average 6.3% statistical and 7.4% systematic 
in the 78 (p, 9) bins - while maximizing the amount of spectral 
information provided by the measurement. It should be noted 
that the magnitude of fractional systematic errors arising from 
the momentum resolution and scale will be affected by the fine- 
ness of the binning. In particular, in the present paper, the mo- 
mentum scale uncertainty has been reduced from 5% to 2% 
since the p-Al publication yet this does not lead to a smaller 
systematic contribution on the measured cross-section. This is 
expected since, simultaneous to the improved reconstruction, 
most momentum bins have been naiTowed by a factor of 2. 

In the end, the statistical plus systematic uncertainty on 
the total integrated cross-section has improved from 5.8% in 
p-Al to 4.9% in p-Be. Due to rebinning and larger statistical 
errors, the average bin-to-bin uncertainty on the differential 
cross-section has changed from 8.2% in p-Al to 9.8% in p-Be. 

We point out that a re-analysis of the proton-aluminum data 
incorporating these changes yields results consistent with those 
published in LIOJ within the systematic errors reported there. 



3 Calculation of the double-differential 
inelastic cross-section 

The goal of this analysis is to measure the inclusive yield of 
positive pions from proton-beryllium collisions at 8.9 GeV/c: 



Be 



X 



The absolutely normalized double-differential cross-section for 
this process can be expressed in bins of pion kinematic vari- 
ables in the laboratory frame, (p^, 6'^), as 
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where: 

2 7r + 

- '^j^^Q is the cross-section in cm^/(GeV/c)/sr for each 
(p^, 6^) bin covered in the analysis 

- j^-^ is the reciprocal of the number density of target nuclei 
for beryllium (1.2349 ■ lO^'^ per cm^). A is the atomic mass 
of beryllium, TVa is Avagadro's number and p is the density 
of beryllium. 

- t is the thickness of the beryllium target along the beam di- 
rection. The target has a cylindrical shape, with a measured 
thickness and diameter of t = (2.046 ± 0.002) cm and d 
~ 3.0 cm, respectively. 

- Ap and Afl are the bin sizes in momentum and solid angle, 
respectively^ 



- A^pot is the number of protons on target after event selection 
cuts (see section |4n i. 

- N'^ {ptt, Ott) is the yield of positive pions in bins of mo- 
mentum and polar angle in the laboratory frame. 

The true pion yield, N^^ {pTr,6T,), is related to the measured 

one, N"^^ {p^, 6*^), by a set of efficiency coiTections and kine- 
matic smearing matrices. In addition, there is a small but non- 
negligible mis-identification of particle types, predominantly 
between pions and protons. Therefore, both yields must be mea- 
sured simultaneously in order to correct for migrations. Eq. [T| 
can be generalized to give the inclusive cross-section for a par- 
ticle of type a 
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Ap — pmax — Pmin; Afl — 2-K{cOs{dr, 
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where reconstructed quantities are marked with a prime and 
, is the inverse of a matrix which fully describes the 

p&ap u a ^ 

migrations between bins of generated and reconstructed quanti- 
ties, namely: laboratory frame momentum, p, laboratory frame 
angle, 9, and particle type, a. In practice, the matrix AI can be 
factorized into a set of individual corrections, as will be done 
here. The reasons for doing this are threefold: 

- Not all efficiencies and migrations are functions of all three 
variables. Particle identification efficiencies and migrations 
do not depend on the angle, 0, and the tracking efficiency 
and momentum resolution are the same for pions and pro- 
tons. 

- Using techniques described below the tracking efficiency 
and particle identification efficiency and migrations can be 
determined from the data themselves and do not rely on 
simulation. This is, of course, preferable wherever possible. 

- Measuring and applying the coiTections separately will ease 
the assessment of systematic eiTors as will be discussed in 
Section|7] 

The form of the corrections can be separated into two basic 
categories: absolute efficiencies and bin-to-bin migrations be- 
tween true and reconstructed quantities. In particular, migra- 
tions in momentum and in particle identification are carefully 
considered. The various efficiency coiTections can, therefore, 
be functions of either reconstructed quantities or true ones, and 
must then be applied at the appropriate point in the analysis. 
This is important given that some corrections, as mentioned 
above, are measured from the data themselves where one has 
only reconstructed quantities. 

Further, we are interested only in secondary 7r+ created in 
primary interactions of beam protons with beryllium nuclei. 
Pions created in interactions other than pH-Be at 8.9 GeV/c are 
a background to the measurement. Tertiary particles are those 
created when secondary particles decay or inelastically interact 
downstream of the target in air or detector materials and are not 
to be included in the measured cross-section. 

In the present analysis, AI~} , „, , has been factorized into 

p6ap 6 a 

the following components. Note that 6^ = ta.n^^ (p^/pz) and 
9y = ta,n~^ {py/pz) are useful variables for viewing the de- 
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lector in x, j/-plane coordinates and are related to the standard Af 



), and the population in each true bin is comprised of 



polar angle by 6 = tan ^(^tan^ 9^ + tan^ 9y). 

- e'''='=°" ,6^,dy) is the efficiency for the reconstruction of 
an 'analysis track'. An 'analysis track' is defined to include 
a momentum measurement as well as a matched time-of- 
flight hit needed for particle identification such that = 

gtrack . gTOFW-match 

- £^^'^{9) is the correction for the geometric acceptance of 
the spectrometer and is a purely analytical function based 
on the assumption of azimuthal symmetry in hadron pro- 
duction and the fiducial cuts used in the analysis. See |10| 
for a full description of the acceptance correction and its 
dependence on the 9y fiducial volume cut. 

- M~]{9 ) is the matrix describing the migration between 

bins of measured and generated momentum. There is a unique 
matrix for each angular bin in the analysis since the mo- 
mentum resolution and bias vary with angle. 

- A/^/ ij)) is a unit matrix, implying that angular migrations, 
which are small, are being neglected. 

- rj'^^^°''^{p,9x,9y,Q) is the absorption plus decay rate of 
secondary particles before reaching the time-of-flight wall 
which is required for particle identification. 

- (1 — 77*°'* (p ,9^,9y,a)) coiTects for the fraction, 77'°'* = 



-tort 

— — , of total tracks passing reconstruction cuts, iV"^ 
which are actually tertiary particles, TV''^'^"*'"^' 



tracks from all reconstructed momentum bins. This approach 
avoids the difficulties associated with inverting a large smear- 
ing matrix due to potential singularity of the matrix as well as 
potential pathologies in the inverted matrix caused by a loss 
of information at the kinematic boundaries of the matrix itself. 
The drawback to this method is that one has some sensitivity 
to the underlying spectrum in the Monte Carlo used to generate 
the matrix (see Sec.|5]l. 

Second, there is a background associated with beam pro- 
tons interacting in materials other than the nuclear target (parts 
of the detector, air, etc.). These events can be subtracted by us- 
ing data collected without the nuclear target in place. We refer 
to this as the 'empty target subtraction': 

N"' {p ,9')^ [7v£,,t {p ,9')- n£^,^ {p,9')]. 



The final form of the cross-section calculation is then given 
by making the above substitution into Eq. [3] 



4 Track selection and reconstruction 
efficiency corrections 



'(p, a) is the efficiency for particles of type a pass- 4.-| Event selection 



ing the electron veto cut used to remove electrons from the 
analysis track sample as described below. 
- M~\ {p) is the particle identification efficiency and migra- 
tion matrix, assumed uniform in 9. 

Once again primed variables are those measured and unprimed 
variables are the true quantities (i.e. after unsmearing). Ex- 
panding Af^^ , into these individual corrections and tak- 

^ ^ pOap 6 a 

ing care of the order in which they are applied gives us the 
final equation for calculating the absolute cross-section from 
the measured yields. 



(Pa" 
dpdQ 
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1 
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NA-p-t ApAQ iVpot 
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, (») : 

0.0. 



X M^^Ap) ■ 



'\p, a) 



1 



1 _ ^absorb 5,^^51^ ^^) (3) 

M-^]{9)-{l~Tf^^\p ,9,,9y,a)) 
1 1 



•7V"'(p,0'). 

There are two additional aspects of the analysis methods 
which are worth mentioning. First, particle distributions are 
built by multiplying a set of correction weights for each recon- 
structed track and weighting events before they are added to the 
total yields. In this way a single reconstructed track is 'spread' 
over multiple true momentum bins according to the elements of 



Protons are identified in the T9 beam at 8.9 GeV/c ex- 
actly as in the 12.9 GeV/c data set and as described in refer- 
ence 1.101 . Two threshold Cherenkov detectors (BCA andBCB) 
placed in the beam line are used to select protons by requir- 
ing a value consistent with the pedestal in both detectors. The 
beam Cherenkov pulse height distributions for the 8.9 GeV/c 
beam are shown in Fig. ID Protons were selected by requiring 
a pulse height less than 120 counts in both detectors, and Fig. 
[3] shows the time-of-flight distributions of those beam tracks 
identified as protons and pions by the Cherenkov selection. 
The beam time-of-flight system is made of two identical scin- 
tillator hodoscopes, TOFA and TOFB, recuperated from the 
previous NA52 experiment and a small target-defining trigger 
counter (TDS). TOFA-TOFB and TOFA-TDS measure time 
differences over a distance of 21.4 m and 24.3 m, respectively. 
We see in Fig.[3]that the two time peaks are consistent with the 
proton and pion hypotheses at 8.9 GeV/c. 

Only events with a single reconstructed beam track in the 
four beam multi-wire proportional chambers (MWPCs) and no 
signal in the beam halo counters are accepted. This MWPC 
track is used to determine the impact position and angle of the 
beam particle on the target. A time measurement in one of three 
beam timing detectors consistent with a beam particle is also 
required for determining the arrival time of the proton at the 
target, <o- This io is necessary for calculating the time-of-flight 
of secondary particles. 

The full set of criteria for selecting beam protons for this 
analysis is as follows: 

- ADC count less than 120 in both beam Cherenkov A and 
beam Cherenkov B 
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Fig. 2. Beam Cherenkov pulse height distributions. BCA in the left panel, BCB in the middle, and BCB vs. BCA in the right panel. The electron 
and pion tagging efficiency is found to be close to 100%; the peaks are separated by ~ 3(7 in both detectors. By requiring a value compatible 
with a pedestal in both Cherenkov detectors the beam protons are clearly separable from pions and electrons as seen in the right panel. 




TOFj - TOF, (ns) TDS - TOP, (ns) 



Fig. 3. Beam time-of-flight distributions. The time difference between 
TOFA and TOFB is shown in the left panel. The right panel is the 
time difference between TOFA and the TDS. The shaded distributions 
are for particles identified as protons by the Cherenkov detectors as 
described in the text. The open histograms are all other beam tracks: 
pions, electrons and muons from pion decays. 



- time measurement(s) in TOFA, TOFB and/or TDS which 
are needed for calculating the arrival time of the beam pro- 
ton at the target, to 

- extrapolated position at the target within a 10 mm radius of 
the center of the target 

- extrapolated angle at the target less than 5 mrad 

- no signal in the beam halo counters 

Prior to the above cuts, for data taken with a nuclear tar- 
get, a downstream trigger in the forward trigger plane (FTP) 
was required to record the eventQ The FTP is a double plane 
of scintillation counters covering the full aperture of the spec- 
trometer magnet except a 60 mm central hole for allowing non- 
interacting beam particles to pass. The efficiency of the FTP is 
measured to be >99.8%. 

Using the FTP as an interaction trigger necessitates an ad- 
ditional set of unbiased, pre-scaled triggers for absolute nor- 
malization of the cross-section. Beam protons in the pre-scale 

^ empty target data sets are recored with an unbiased trigger setting 
since these samples are used to calibrate the experimental apparatus 
and not just in the empty target subtraction for cross-section measure- 
ments. 



trigger sample (1/64 of the total trigger rate for the 8.9 GeV/c 
Be data set) are subject to exactly the same selection criteria 
as FTP trigger events allowing the efficiencies of the selections 
to cancel and adding no additional systematic uncertainty to 
the absolute normalization of the result. These unbiased events 
are used to determine the A'pot used in the cross-section for- 
mula and listed in Table 14.21 The number of protons-on-target 
is known to better than 1%. 

Applying these criteria we are left with the event totals 
summarized in Table 14.2] 



4.2 Secondary track selection 

The following criteria have been applied to select tracks in 
the forward spectrometer for the accepted events: 

- a successful momentum reconstruction using downstream 
track segments in NDC modules 2, 3, 4 or 5 and the posi- 
tion of the beam particle at the target as an upstream con- 
straint (here, upstream and (iownstream are relative to the 
spectrometer magnet); 

- a reconstructed vertex radius (i.e. the distance of the recon- 
structed track from the z-axis in a plane perpendicular to 
this axis at z = 0) r < 200 mm; 

- number of hits in the road around the track in NDCl > 
4 and average for these hits with respect to the track in 
NDCl < 30 (this is applied to reduce non-target interaction 
backgrounds); 

- number of hits in the road around the track in NDC2 > 6 
(this is applied to reduce non-target interaction backgrounds); 

- a matched TOFW hit passing the quality cuts described in 
Sec. 11231 

- reconstructed angles are within the fiducial volume to be 
used for this analysis, —210 mrad < 9x < mrad and - 
80 mi-ad < 61^ < 80 mrad. 

These cuts are identical to those used in the analysis of the p-Al 
data except for the reconstructed vertex radius < 200 mm cut. 
It was found that due to a feature of the algorithm this addi- 
tional requirement improved the momentum resolution consid- 
erably at reconstructed momenta below w 1 .5 GeV / c. 
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Data Set 


Be 5% 8.9 GeV/c 


8.9 GeV/c Empty Target 


protons on target 


13,074,880 


1,990,400 


total events processed 


4,682,911 


413,095 


events with accepted beam proton 


2,277,657 


200,310 


beam proton events with FTP trigger 


1,518,683 


91,690 


total good tracks in fiducial volume 


95,897 


3,110 



Table 1. Total number of events in the 8.9 GeV/c beryllium 5% Ai target and empty target data sets, and the number of protons on target as 
calculated from the prescaled trigger count. 



a 1.2- 



• Data 

-MC 



-21 0< 
-80 < ( 



< mrad 
80 mrad 



6 7 8 

p (GeV/c) 



a 1.2 
to 



' Data 

-MC 



0.75 < p <6.50G6V/c 
-80 < e. < 80 mrad 



-0.2-0.15-0.1-0.05 -0 0.05 0.1 0.15 0.2 

8, (rad) 





• Data 




MC 




: 0.75 < p < 6.50 GeV/c 




: -210 < < mrad 





, (rad) 



Fig. 4. Tracking efficiency for positive particles traversing the detector using the target as upstream track constraint as a function of particle 
momentum (upper left), production angle in the horizontal plane, 0x (upper right), and production angle in the vertical plane, 9y (lower). The 
9y plane is orthogonal to the spectrometer bending plane and not sensitive to the momentum dependent acceptance, and the bottom panel shows 
the purest measure of the average track reconstruction efficiency within the fiducial volume to be 96%-97%. 



Applying these cuts to reconstructed tracks in accepted even 
we are left with 95,897 total good tracks in the beryllium thin 
target data set as listed in Table |4.2] 



4.3 Track reconstruction efficiency 



The track reconstruction efficiency has been measured from 
the data exactly as described in ifTOl . The efficiency is shown 
in Fig. m The effects of the two changes in track reconstruc- 
tion are evident in the efficieny curves. First, the efficiency is 
now flat and ~ 97% above 2 GeV/c due to the improvement 
in the minimization done as part of the tracking algorithm. 
Second, the loss of efficiency at momenta below 1 .5 GeV /c is 
due to the reconstructed vertex radius < 200 mm cut discussed 
above. But, as before, the tracking efficiency can be measured 
from the data themselves, so the systematic error on the correc- 
tion comes only from the statistical uncertainty in the sample 
used to calculate the correction. As in the p-Al publication, the 
12.9 GeV/c aluminum data and the 8.9 GeV/c beryllium data 
have been combined to minimize this uncertainty. The drop in 
efficiency at large, positive values of is due to geometric ac- 
ceptance as low momentum tracks are bent out of the spectrom- 
eter missing the downstream chambers. The present analysis is 
performed using tracks in the range —0.210 rad < 0^ < rad 
where the acceptance is flat in momentum. 



4.4 Absorption, decay and tertiary track corrections 

The correction for absorption and decay refers to secondary 
particles created in the nuclear target that never make it to 
the time-of-flight wall for detection and possible identification. 
Figure [T] shows the location of the time-of-flight scintillator 
wall just beyond the back plane of drift chamber modules NDC3, 
NDC4 and NDC5. We use the Monte Carlo simulation to de- 
termine the size of the correction and the result is shown in 
Fig- 13 Note this is an upward adjustment to the raw yield mea- 
sured and is implemented as 1/(1 — ri"'^''°^^{p,9x,0y,a)) in 
Eq.[3] The absorption correction (which includes pion decays) 
is a function of 9^ and 9y because it depends on the amount 
and type of physical material a particle passes through, thus 
the geometry of the detector. It is a function of a because of 
the different interaction cross-sections and possible decay rates 
of hadrons. This correction is separated from the tertiary cor- 
rection discussed below because it does not depend on event 
multiplicity, kinematics or other details of the hadron produc- 
tion model used in the simulation, but only the total interaction 
cross-sections which are significantly more certain. In fact, the 
relevant cross-sections are typically known to « 10% and we 
assume this uncertainty on the absorption correction just as in 
the p-Al publication. Because the correction is of order 30- 
40% for pions, the average systematic error contribution to the 
cross-section turns out to be 3.6%. 

The tertiary correction refers to the subtraction of recon- 
structed tracks which are actually reconstructions of tertiary 
particles, i.e. particles produced in inelastic interactions or de- 
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Fig. 5. Absorption corrections for pions and protons according to Monte Carlo simulation as a function of particle momentum (left), production 
angle in the horizontal plane, 6x (center) and production angle in the vertical plane, 6y (right). 
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Fig. 6. Tertiary particle rates for pions, muons (which get identified as pions) and protons according to Monte Carlo simulation as a function 
of particle momentum (left), production angle in the horizontal plane, 0x (center) and production angle in the vertical plane, 9y (right). 



cays of true secondary particles and not in primary interac- 
tions of 8.9 GeV/c protons with beryllium nuclei (See sec- 
tion O. The tertiary subtraction includes muons created in de- 
cays which are falsely identified as pions nearly 100% of the 
time due to their high (3. The correction is significantly smaller 
than the absorption correction (compare Figs. |5] and |6]l, but is 
less certain, so the contribution to the systematic eiTor is non- 
negligible. This tertiary subtraction is also generated using the 
Monte Carlo simulation but is dependent on the details of the 
hadron production model used in the simulation. Most of the 
material where tertiary particles might be produced in the de- 
tector is carbon, so it is the simulation of inelastic interactions 
of low-energy protons and pions in carbon that become im- 
portant in generating this correction. Previously this correction 
was assumed to be 100% uncertain, but comparisons of low 
momentum HARP p+C, tt^+C and tt^+C data to the hadronic 
models used in the simulation have verified these models to 
« 50% and allowed us to lower the systematic error on this 
correction. Fig. |6] shows the average size of the correction to 
the 7r+ yield to be about 4% (2% 7r+ + 2% so the average 
systematic error on the cross-section coming from the subtrac- 
tion of tertiaries ends up being 1.8%. 



5 Momentum resolution and scale 
corrections 

Two important sources of systematic uncertainty are the 
momentum resolution and the absolute momentum scale of the 
reconstruction. Each, in general, vary with the value of the true 
momentum and angle of the track. The Monte Carlo will be 
used to generate corrections to the measured momenta so a 
validation of the simulation becomes necessary. Three tech- 
niques have been developed to compare the resolutions of re- 
constructed quantities in data and Monte Carlo of which only 
one was available at the time of the p-Al paper. The challenge 
is to isolate a set of tracks in the data sample with a known 
momentum. The three methods are described fully below but, 
briefly, they are based on empty target data sets, samples of 
elastic scattering events and using the excellent resolution of 
the time-of-flight system to determine the momentum. Results 
from these studies have been used to tune the detector simula- 
tion and to estimate a systematic uncertainty due to momentum 
resolution and scale. 

A multiplicative momentum scale coiTection is applied to 
all reconstructed tracks in the data to remove a 6^, Oy depen- 
dence seen in calibration samples. After this correction we see 
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no significant momentum mis-calibration beyond the 2% ab- 
solute momentum scale uncertainty estimated using the elastic 
scattering technique (see Sec. l5.2l l. The minimum ratio of mo- 
mentum bin width over momentum bin central value is 8%, 
four times this value. An average uncertainty on the 7r+ cross- 
section of 3.6% is estimated from this effect (the effect is de- 
pendent on the steepness of the pion and proton spectra, the 
size of the momentum bias and the momentum binning). 

The resolution of the momentum measurement leads to smear- 
ing in the reconstructed pion spectrum. Momentum 'unsmear- 
ing' is performed in the analysis by using the Monte Carlo 
to generate a momentum migration matrix, M^]{0 ) which 
describes how a reconstructed momentum value is distributed 
across bins of true momentum. The accuracy of this migration 
matrix, which is generated from Monte Carlo samples, depends 
upon two factors - the similarity of the underlying hadron dis- 
tributions to those in the data and the agreement between the 
resolutions of reconstructed quantities in the data and Monte 
Carlo relative to the true momenta. The latter, the agreement 
between the resolutions of reconstructed quantities in data and 
Monte Carlo, is the topic of the next three subsections. The 
techniques developed have been used to tune parameters of the 
simulation to achieve good agreement. The first effect, the sim- 
ilarity of the hadron spectra to the data, therefore, will domi- 
nate the systematic error for this correction. We estimate a mo- 
mentum resolution uncertainty by performing the cross-section 
analysis using migration matrices generated from Monte Carlo 
samples which use different hadronic models. The average un- 
certainty on the 7r+ cross-section estimated from this effect is 
3.4%. 

5.1 Momentum calibration using empty target data 
sets 

The first method for calibrating the momentum reconstruc- 
tion uses empty target data sets where the incoming beam mo- 
mentum value acts as the known momentum. Fig.|2]shows the 
result of such a study using empty target data and Monte Carlo 
samples for beam momenta of 1.5, 3.0, 5.0, 8.0, 8.9, 12.0, 12.9 
and 15.0 GeV/c. The shape of the (Jp/p curve versus momen- 
tum is as expected, and the agreement between data and Monte 
Carlo is excellent. For the previous analysis |10| only these 
'test-beam data' were available, sampling the spectrometer re- 
sponse on the z-axis only. Now, additional methods have been 
used. A method using samples of elastic scattering events ex- 
tends the range to larger angles while a method employing the 
time-of-flight detector covers the region of low momenta. Thus 
the full range of relevant angles and momenta is covered. 

5.2 lUlomentum calibration using elastic scattering 
events 

The elastic scattering process provides one track in the for- 
ward direction with a momentum close to the beam momen- 
tum and a soft large angle proton. In p-p and 7r-p scattering 
the kinematics is fully determined by the measurement of the 
direction of one of the outgoing particles. The precision of the 
measurement of the angle of the forward scattered particle is 
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Fig. 7. Momentum resolution a{p)/p (Gaussian fit) as a function of 
momentum (in GeV/c) for the drift chambers: the data were taken 
using several well-defined discrete beam momenta and no target for 
pions (open circles) and protons (open squares). Also shown (filled 
circles) is the corresponding resolution found using the Monte Carlo 
simulation. 



sufficient to predict the momentum of that particle without sig- 
nificant error. Elastic events can be readily selected by impos- 
ing combined criteria in the large angle and forward spectrom- 
eter The main selection is the requirement of one and only one 
track in the TPC, identified as a proton, and exactly one track in 
the forward direction. Further constraints have been used, such 
as a match of the kinematics of these tracks and the selection 
of events with exactly the expected number of hits in the trig- 
ger counters and the barrel RPC detectors around the TPC. In 
these selections no constraints on the momentum measurement 
of the forward track has been set. The purity of the selection 
of the elastic scattering events used in this analysis can be es- 
timated by studying the p-9 distribution of the proton recoil 
tracks measured in the TPC. The two quantities are fully cor- 
related for elastic scattering events. From the small number of 
events outside the expected region, one can estimate the back- 
ground in the sample not to exceed 1%. The systematic error 
introduced by this background is negligible compared to the 
2% estimated overall momentum calibration error. 

Figure [8] shows the results of an analysis of elastic events 
using 3 GeV/c, 5 GeV/c, and 8 GeV/c beams impinging on a 
hydrogen target. Incoming pion and proton data are combined. 
The left panel reveals a momentum offset in the data of about 
2% at all momenta and in three angular regions, 30 mrad- 
60 mrad, 60 mrad-100 mrad and 100 mrad-150 mrad. This 
2% will contribute a systematic uncertainty to the final cross- 
section and will be discussed in Sec. 17. II The right panel shows 
the RMS of the momentum measurement as a function of mo- 
mentum in the same two angular regions. The resolution mea- 
sured from data is compared to that in the Monte Carlo and, as 
with the empty target samples, the agreement is excellent. 

The momentum scale and resolution measured with the elas- 
tic events extends the calibration toward larger angles than probed 
with the empty target data alone and make it possible to char- 
acterize the spectrometer over a larger range of its aperture. 
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Fig. 8. Momentum scale A{p) /p and resolution cr(p) /p (Gaussian fit) 
as a function of momentum in three angular regions (circles: 30 mrad- 
60 mrad; boxes: 60 mrad-100 mrad; diamonds 100 mrad-150 mrad). 
The data were taken using 3 GeV/c, 5 GeV/c and 8 GeV/c beams 
impinging on a hydrogen target. The top panel shows no significant 
momentum mis-calibration beyond the 2% absolute momentum scale 
uncertainty estimated using the elastic scattering technique. The bot- 
tom plot shows the resolutions from data (open symbols) compared 
with the corresponding Monte Carlo simulations (filled symbols). In- 
coming pion and proton data are combined. 



However, the lowest momentum at which these data are avail- 
able is 3 GeV/c, such that another method had to be developed 
to study the performance of the spectrometer for lower momen- 
tum particles. 



5.3 Momentum calibration using time-of-flight 

The TOFW system fT2| can be used to provide a momen- 
tum calibration in a momentum range where the dependence of 
the time-of-flight on the momentum for pions is much smaller 
than for protons and at the same time the resolution of the 
TOFW is better than the prediction of f3 for protons based 
on the momentum measurement. For protons the sensitivity 
to the momentum resolution is larger than all other dependen- 
cies in the beta resolution in the region 0.22 GeV/c < p < 
2.0 GeV/c. The TOFW/3 resolution is typically cr(/3)//3 = 
0.005. The precision o'(/3pi-ed)//^ ™ the calculation of /3 for 
protons ranges from 0.05 at 0.22 GeV/c to 0.005 at 2 GeV/c 
due to the momentum resolution which is of the order of 
100 MeV/c in this range. Therefore the width of the /? peak for 
a sample of protons selected in a small range in measured mo- 
mentum shows a large sensitivity to the momentum resolution. 

To exploit this feature for the determination of the mo- 
mentum resolution at small momenta the time-response of the 
TOFW has to be measured for pions and protons separately. 
First a very clean sample of pions is selected. Particles of nega- 
tive charge are selected for this purpose to provide negative pi- 
ons. In principle, this sample can be contaminated by electrons 
and negative kaons. Antiprotons are expected to be negligible. 
At a momentum below the Cherenkov threshold for pions, elec- 
trons are rejected by retaining only particles without a signal in 
the Cherenkov detector. The remaining sample contains mainly 
TT^ with a small background of K^. This kaon background is 
visible in the TOF spectrum and has been taken into account. 




^ * 



Fig. 9. RMS momentum resolution as a function of momentum in four 
angular bins (circles: 30 mrad-60 mrad; boxes: 60 mrad-100 mrad; 
diamonds: 100 mrad-150 mrad; triangles: 150 mrad-200 mrad). Data 
are shown with open symbols; Monte Carlo with filled symbols. 



At a momentum above the pion Cherenkov threshold, the 
TT^ sample is selected by requiring a signal in the Cherenkov. 
Although the electron background is expected to be small at 
these momenta, a selection with the calorimeter is used to re- 
ject them. The same selection is also used at momenta be- 
low Cherenkov threshold in addition to the rejection with the 
Cherenkov detector. Thus a sufficiently clean sample of nega- 
tive pions can be obtained in the whole momentum range. 

From the sample of positive particles, positrons are removed 
by the selection using the electron identifier Below Cherenkov 
threshold only particles without signal in the Cherenkov are 
used. This selection retains protons, kaons and pions below 
pion Cherenkov threshold and only protons and kaons above 
the threshold. The precision of the TOFW is sufficient to pro- 
vide a good separation of pions and protons below pion 
Cherenkov threshold. Thus, in the whole range relevant for the 
experiment a clean sample of protons can be obtained, albeit 
with a small contamination of kaons. As mentioned above, the 
measurement with the negative pions characterizes the TOFW 
response, both its absolute time and its time resolution. The 
negative pions also provide a perfect prediction for the behav- 
ior of the TOFW measurement for the positive pions. 

The measurement of the properties of the momentum de- 
termination in the spectrometer is then obtained by selecting 
small regions (bins) of measured momentum and fitting the (3 
spectrum of protons with a function which takes into account 
the width of the momentum bin, the calibrated f3 resolution and 
as a free parameter the momentum resolution. The results for 
data and Monte Carlo and in four angular regions are shown 
in Fig. |9] The resolution measured should be interpreted as an 
RMS of the momentum resolution and is larger than the cr of 
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the Gaussians fitted to the direct beam data using the runs with- 
out target, but consistent with the RMS of the latter. 



Combining information from these three techniques one is 
able to map out the momentum resolution and scale in both data 
and Monte Carlo for comparison. The results indicate good 
agreement across a range of momenta and angles allowing us to 
utilize the Monte Carlo simulation to generate the momentum 
resolution correction matrix, M~] {9 ). 



6 Particle identification 

The particle identification method used here follows ex- 
actly that used in the analysis of the p-Al data 1 10|. However, 
a significant improvement in efficiency and purity of the PID 
and, consequently, a large reduction in the systematic uncer- 
tainties on the cross-section have been realized by improving 
the association of PID detector hits with reconstructed tracks. 
We include here a description of the PID method for complete- 
ness and clarity and a full explanation of how detector hits are 
selected follows in the next section. 



6.1 The pion-proton PID estimators and PID 
efficiency calculation 

The analysis uses particle identification information from 
the time-of-flight and Cherenkov PID systems; the discrimina- 
tion power of time-of-flight below 3 Ge V / c and the Cherenkov 
detector above 3 GeV /c are combined to provide powerful sep- 
aration of pions and protons. The calorimeter is presently used 
only for separating pions and electrons when characterizing the 
response of the other detectors. The resulting efficiency and pu- 
rity of pion identification in the analysis region is excellent. 

Particle identification is performed by determining the prob- 
ability that a given track is a pion or a proton based on the ex- 
pected response of the detectors to each particle type and the 
measured response for the track. Information from both detec- 
tors is combined for maximum discrimination power using a 
Bayesian technique. 



P(a|/3,7Vpe,p,0) = 



Pil3,Np,\a,p,9)-Pia\p,0) 
E PW,Np,\i,p,e)-P{i\p,9) 

-TT.p.... 



(4) 



spectra; that is, the prior distributions will be flat and equal ev- 
erywhere, P{i\p,9) — 1 for all p, 9. This allows the priors 
to be dropped from the expression, but the PID estimator no 
longer has a full probabilistic interpretation and cannot be di- 
rectly used to estimate the particle yields. One could iterate the 
probability distributions to determine the yields. Alternatively, 
one can build the PID estimator for each track independently, 
and an efficiency and migration must be determined for a given 
cut on the estimator value, Ptrack > ^cut- We will see that the 
necessary corrections are small and the systematic uncertainty 
is negligible compared to other sources, making this approach 
adequate. Second, we will consider the response functions of 
the different PID detectors as independent and can therefore 
factorize the probability into separate terms for the TOFW and 
CHE. Third, with the new detector hit selections, the time-of- 
flight and Cherenkov detector responses show no angular de- 
pendence allowing 9 to be removed from the above expression. 
Finally, we will only consider pions and protons as possible 
secondary particle types. Monte Carlo simulation shows other 
potential backgrounds to a 7r+ yield measurement to be small. 



P(a|/3,iVpe,p) 



P(/3|a,p)-P(iVpe|a,p) 



(5) 



P(/3, |7r,p) • P(7Vpo|7r,p) + P(/3|p,p) • P(iVpe|p,p) 



where P{a\i3, Np^p) is the PID estimator for a track with re- 
constructed (3, Npc and p to be of type a and P(/3|p,p) and 
P(iVpo|p, p) are the response functions for the TOFW and CHE, 
respectively, which will be fully described below. Pions are se- 
lected by making a cut in this PID variable, equal to 0.6 in the 
present analysis. 

The efficiency for pion selection and the migration between 
pions and protons can be calculated analytically from the parametrized 
detector response functions for a given cut in probability. These 
corrections consist (for each momentum bin) of a 2x2 matrix 
(PID efficiency matrix) 



M^TT M^p 

AdpTi A/pp 



(6) 



where P(a|/?, Npc, Pi is the probability that a track with re- 
constructed velocity /3, number of associated photo-electrons 

Npc, and momentum and angle p and is a particle of type 7\f — C^' ■ I 

a. P{i\p,9) is the so-called prior probability for each parti- g'^^ i -Jb 

cle type, i, and is a function of p and 6. In the Bayesian ap- 
proach, the priors represent one's knowledge of the relative 
particle populations before performing a measurement. Finally, 
P(/3, Npc\i,p, 9) is the expected response (/3 and A^pc) of the 
PID detectors for a particle of type i and momentum and angle 

The following simplifications are applied to Eq.|4] First, we 
will assume no a priori knowledge of the underlying pion/proton 



where Mtttt and Mpp are the efficiencies for correctly identi- 
fying pions and protons, respectively, and M-^p and Mp-n are 
the migration terms of true protons identified as pions and true 
pions identified as protons, respectively. The full covariance 
matrix (16 terms) of the PID efficiency matrix is also com- 
puted analytically. This calculation is explained fully in ifTTl . 
The final expression for the PID efficiency-migration matrix 
elements reads: 



d/3 G,(/?) + ■ S{Pf^ IcHE > Pent) 



(7) 



where Cj" is the Cherenkov efficiency (Fig. [15]). Gj{l3) and 
wj^ are the Gaussian (Fig.[T2]i and non-Gaussian contributions 

to the beta response, respectively, and S{Pf "{cue > Pcut) is 
a step function controlling the integration of the non-Gaussian 
part, uj. 
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The 2x2 matrix of Eq. |6] is easily inverted for converting 
reconstructed yields into true yields of pions and protons. The 
elements of the matrix in Eq.|6]are shown in Fig.[T6]as a func- 
tion of particle momentum. The pion efficiencies are >95% 
and the proton-pion migrations are all less than 1%. But first 
we must desribe the detector response functions that are a key 
input to the efficiency and migration calculations. 

6.2 PID detector hit selection and response 
functions 

This section describes the quality criteria applied to select 
PID detector hits and the resulting response functions for pi- 
ons and protons. Previously, PID detector hits have been asso- 
ciated with reconstructed tracks based only on a geometrical 
matching criterion. A Kalman filter package is used to extrap- 
olate each track to the plane of each detector, and this posi- 
tion is compared to the reconstructed x, y positions of all re- 
constructed hits in that detector The detector hit with the best 
matching was then assigned to that track. In this scheme the 
reconstruction can associate a single detector hit with multiple 
tracks and each track likely has additional candidate detector 
hits which are being ignored. In particular, in |10| it was seen 
that a fraction of protons had a non-negligible amount of as- 
sociated photo-electrons due to light from pions or electrons 
being wrongly associated with proton tracks. Also in ifTOI the 
TOFW response function contained a non-gaussian component 
where «10% of reconstructed (3 fell greater than Scr from the 
mean expectation and could not be used for identification. Ad- 
ditional criteria have been developed for selecting PID detector 
hits to address these issues and have led to a reduction in PID 
backgrounds by as much as a factor of 10 in some regions of 
phase space. This is the source of the drastic improvement in 
PID systematics since our previous publication (3.5% to 0.5%). 
The PID hit selection is described in detail below. 



6.2.1 Time-of-flight response 

A time-of-flight measurement is required for particle iden- 
tification in this analysis. It was discovered, due to the presence 
of a significant, almost flat background far from the Gaussian 
peaks in the /3 distributions, that a more strict set of selection 
cuts was required to ensure a quality time-of-flight measure- 
ment. The small efficiency loss due to this selection can be 
measured directly from the data and will be combined with the 
tracking efficiency discussed above to form an overall recon- 
struction efficiency, 

^rccon ^track _ ^TOFW — match 

Each track can have multiple time-of-flight measurements 
(TOFW — to) associated with it in the reconstruction. It is pos- 
sible for a single hit to match with multiple tracks if the tracks 
are close enough together when hitting the wall. It is also pos- 
sible that electromagnetic showers associated with a particle 
passing through detector material can create additional hits be- 
yond the primary hit caused by the hadron of interest. To min- 
imize inaccurate time measurements due to these effects, the 




r of hit matching mips (PMT1+PMT2) 



Fig. 10. TOFW hit reconstructed variables. The left panel shows the 
distribution of the between the extrapolated track position and the 
reconstructed scintillator hit position. The right panel shows the total 
reconstructed number of minimum ionizing particles (mips) from the 
two PMTs on the scintillator volume that was hit. Time-of-flight hits 
are selected by requiring a x'^ < 6 and number of mips > 1.5 (see the 
text). 



time-of-flight candidates for each track are time-ordered and 
the earliest hit passing the following criteria is selected: 

- if the track shares the TOFW hit with another track, it must 
have the better geometric matching ; 

- of the geometrical matching between the track and TOFW 
hit < 6; 

- total reconstructed number of minimum ionizing particles 
(mips) from the two PMTs in a hit > 1.5. 

The distribution for track-TOFW hit matching and the total 
pulse-height distribution for this data set are shown in Fig.fTol 

Having applied these criteria to time-of-flight measurements 
we must understand the associated efficiency loss as well as the 
remaining level of non-Gaussian component to the f3 spectrum. 
Each of these have been carefully measured and the needed 
corrections applied. 

The efficiency is measured from the data by using a sample 
of reconstructed tracks which leave a signal in the calorime- 
ter (downstream of the scintillator wall) and asking how often 
a time-of-flight measurement passing selection cuts is found. 
Fig.[TT|shows the matching efficiency for data and Monte Carlo 
to be flat in momentum and around 95% in the data. This effi- 
ciency is combined with the tracking efficiency to provide the 
total analysis track reconstruction efficiency. The efficiencies 
generated from the data themselves have been used in the anal- 
ysis for the results being presented here. 

Having applied these selection criteria we must character- 
ize the response of the time-of-flight detector for different par- 
ticle types. The (3 response has been parameterized by a Gaus- 
sian function. A method has been developed to extract the pa- 
rameters of this function directly from the data and has aheady 
been described in |10| and |11|. The method was applied to 
both data and Monte Carlo events and the results are shown 
in Fig. [12] The figure shows a small bias in the simulation of 
the time-of-flight resolution for protons above « 2 GeV /c, re- 
sulting in a particle identification efficiency bias of about 0.5% 
(see Fig.[T6]l. To avoid any bias introduced by the simulation of 
the time-of-flight system, the response functions as determined 
from data are used in the analysis of data. 
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Fig. 11. TOFW matching efficiency as a function of particle momen- 
tum (left) and production angle in the horizontal plane, 9x (right) as 
measured from data and Monte Carlo. The TOFW inatching efficiency 
does not have the momentum dependence of the tracking efficiency, 
but does exhibit the same effects of geometric acceptance as the drift 
chambers as seen in the right, 9x, plot. Note the present analysis is 
performed using tracks in the range —0.210 rad < 9^ < rad where 
the acceptance is flat in momentum. 



ttti 



Fig. 13. Pion (left) and proton (right) /3-outlier rates as a function of 
momentum. The pion outlier rate has been measured from the data; 
the proton outlier rate is estimated using the Monte Carlo. In refer- 
ence 111] it has been demonstrated that no bias is seen when using the 
Monte Carlo to calculate the outlier rate. 



6.2.2 Cherenkov response 



The Cherenkov detector is used to veto electrons below 3 
Additionally, there is a small rate of non-Gaussian time- GeV/c and to differentiate pions from protons above 3 GeV/c. 
of-flight measurements, called "/3-outliers", which must be ac- Below 3 GeV/c the Cherenkov signal is not used in the 
counted for separately. /3-outliers are defined as time-of-flight calculation of the particle identification probabiUty according 
measurements greater than 5a from the mean of the expected f3 to Eq. E but instead electrons are removed by a simple veto 
response function. This small, non-Gaussian component of the of tracks with greater than 15 photo-electrons. Fig. [Hdemon- 
time-of-flight response, shown in Fig. [13] has been fully ac- strates the effect of the electron veto cut. The left panel shows 
counted for in the PID efficiency calculation described above, the e/n ratio in the Monte Carlo before and after applying the 
It should be noted that, due to the improvements in the hit se- 15 photo-electron cut below 3 GeV/c. The remaining elec- 
lection criteria being described here, the /3-outlier effect as de- tron contamination is less than 1% everywhere, and less than 
scribed in 1 10| has been reduced from -10% to -1-3% since 0.5% in the region where the veto is applied. One expects a 
that publication. The outlier rate is the largest contribution to very small efficiency loss for pions and protons due to this cut 
the systematic error coming from particle identification, and in photo-electrons and this is also shown in Fig. [H Approxi- 
with this improvement PID now makes a negligible contribu- mately 1% of pions and protons do not pass the electron veto; 
tion to the total systematic error in the cross-section analysis, a correction has been applied in the present analysis. 

Above 3 GeV /c the Cherenkov is a powerful discriminator 
of pions and protons. (Monte Carlo simulations indicate that 
there are a negligible number of electrons above 3 GeV/c and 
these are thus ignored.) Presently the Cherenkov is being used 
digitally. That is the spectral information of the light output is 
not being used. Instead we define a signal as an associated hit 
protons with greater than 2 photo-electrons. Two or less is considered 

~ no signal. Based on this definition we determine the efficiency 

1— for pions and protons to have a signal in the Cherenkov as a 

■ — function of particle momentum. Fig.[T5]shows the expected re- 

ooDs- pionr^"~~"'™™"~^-^ sponse for pions and protons in the Cherenkov both above and 

below the pion threshold. Above threshold the Cherenkov is 
05 1 1.5 2 2.5 3 greater than 99% efficient for pions. The small efficiency for 

p(GeV/c) p(GeV/c) ^ r J 

protons and pions below threshold of around 1.5% is due to 

TP- It ^ f ■ jfjf false associations with light generated by Other particles in the 

Fig. 12. p response for pions and protons as measured from data fe fe j f 

and Monte Carlo. The left panel shows the mean beta values with event. 

the error bars representing the width of the Gaussian distribution. The 

right panel highlights the width of the distributions and shows how 

the resolution asymptotically approaches ^ 0.006. The solid points Using the characterized responses of the TOFW and CHE 

are the response as measured from data; the dashed histograms are the detectors we can calculate the PID estimator for reconstructed 

response as measured from Monte Carlo. tracks given in Eq. |5] and the efficiency-migration matrix ele- 
ments given by Eq.|7]and shown in Fig.fTSI 
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Fig. 14. The left panel shows the e/iv ratio from a Monte Carlo simu- 
lation before (solid points) and after (open squares) the application of 
a 15 photo-electron cut. This cut reduces the electron contamination 
to 0.5% or less in the region where it is applied. The right panel shows 
the efficiency for pions (solid points) and protons (open squares) to 
pass the 15 photo-electron cut below 3 GeV/c, and is ~ 99% for 
both. 





Fig. 16. Particle ID efficiency and migration matrix elements as a 
function of momentum. Upper left is the pion identification efficiency, 
Af^^. Upper right is the proton to pion migration, A/,rp. Lower left is 
the pion to proton migration, A/p^. Lower right is the proton iden- 
tification efficiency, A/pp. Points with errors are the values calcu- 
lated from data. The dashed histograms are those determined from the 
Monte Carlo. The slight data-Monte Carlo bias seen above 2.5 GeV/c 
results from the bias seen in the TOFW simulation. The values deter- 
mined from data have been used in the analysis to avoid sensitivity to 
this bias. 
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Fig. 15. Cherenkov response for pions (left) and protons (right). The 
points are the efficiencies for a track to have an associated Cherenkov 
hit with greater than 2 photo-electrons. The threshold for pions at 
around 2.6 GeV/c is clearly visible (note the log scale). The small 
efficiency for protons and below threshold pions of around 1.5% is 
due to false associations with light generated by other particles in the 
event. 



7 Physics results 



Applying corrections to the raw yields as described in the 
previous sections and according to Eq. [3] we have calculated 
the double-differential inelastic cross-section for the produc- 
tion of positive pions from collisions of 8.9 GeV/c protons 
with beryllium in the kinematic range from 0.75 GeV/c < 

< 6.5 GeV/c and 0.030 rad < 6^ < 0.210 rad. 

Systematic errors have been estimated and will be described 
below. A full (13 x 6)'^ = 6048 element covariance matrix 
has been generated to describe the coiTelation among bins. The 
data are presented graphically as a function of momentum in 30 
mrad angle bins in Fig.[T7]and ID projections onto the momen- 
tum and angle axes are shown in Fig. [18] The central values and 
square-root of the diagonal elements of the covariance matrix 
are listed in Table El 



7.1 Error estimation 

A full systematic eiTor evaluation has been performed on 
these data in order to estimate the accuracy of the measure- 
ment being presented. Statistical errors from both the beryl- 
lium target data and the 8.9 GeV/c empty target data set used 
to subtract non-target backgrounds are also included. 

The uncertainties associated with the various corrections 
applied have been estimated through a combination of analyt- 
ical and Monte Carlo techniques. The approach used here has 
largely followed the methods of 1 10|. There are eight sources of 
systematic uncertainty considered which can be grouped into 
three basic categories: track yield corrections, PID, and mo- 
mentum reconstruction. 

Only the PID uncertainties have been calculated analyti- 
cally. The covariance matrices of the PID efficiency-migration 
matrices described in Sec.|6]have been calculated and the er- 
rors propagated. The cross-section uncertainties from the other 
sources are estimated by performing the cross-section calcula- 
tion N times for N variations of each correction applied. The 
fully coiTelated error matrix for each correction is then built 
from the N cross-section results. 
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Fig. 17. Measurement of the double-differential production cross- 
section of positive pions, d^a^ /dpdfi, from 8.9 GeV/c protons on 
beryllium as a function of pion momentum, p, in bins of pion angle, 
6, in the laboratory frame. The error bars shown include statistical er- 
rors and all (diagonal) systematic errors. The dotted histograms show 
the extended Sanford-Wang parametrization of Eq.[T4]with parameter 
values given in Table[3] 
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Fig. 18. Projections of the double-differential cross-section results 
onto the momentum axis integrated over the angular range 30 mrad < 
9 < 210 mrad (left) and onto the angle axis integrated over the mo- 
mentum range 0.75 GeV/c < p < 6.5 GeV/c (right). Projections for 
the best-fit extended Sanford-Wang parametrization are also shown, 
as indicated by dotted histograms. 



where / and j label bins of (p, 6), Ef^ is the element of 
one of the error matrices (labeled a), (P a'^y / {dpdf2) is the 
central value for the double-differential cross-section measure- 
ment and cf(J^ ^J{dpdf]) is the cross-section result from the 
n"' variation for the a*'* systematic. For example, to estimate 
the cross-section uncertainty arising from the absorption cor- 
rection, 100 analyses are performed where only the absorption 
correction is randomly fluctuated 100 times with an RMS of 
10%. This 10% is the uncertainty of the absorption correction 
as described in Sec. 14.41 The total error matrix is just the sum 
of the a matrices, Eij — -^fj- 

The full 78 x 78 elements of the covariance matrix will 
not be published here, but to characterize the uncertainties on 
this measurement we show the square-root of the diagonal ele- 
ments of the covariance matrix plotted on the data points in Fig. 
[Tt] Fig. [19] shows the fractional uncertainty (again diagonal) 
for each (p, 6) bin. The total error as well as the contributions 
from statistical errors, track yield corrections and momentum 
reconstruction are shown. Additionally, we define a dimension- 
less quantity, Jdiff, expressing the typical diagonal error on the 
double-differential cross-section 



^diff = 



(10) 



We also define the fractional error on the total integrated pion 
cross-section in the range of the measurement (0.75 GeV/c < 
p < 6.5 GeV/c, 30 mrad < 6* < 210 mrad), Sint- 



5\nt = 



{dpdQ)iEij {dpdQ)j 



(11) 



where {d'^a^)i is the double-differential cross-section in bin i, 
{d^a'^ / {dpdfi))i, multiplied by its corresponding phase space 
element {dpdQ)i. Eij is the covariance matrix evaluated for 
the double-differential cross-section data. 

For example, the 30%-40% absorption correction (see Fig. 
|5]l with a 10% uncertainty results in an average diagonal error 
of 3.6% on the cross-section, as one expects. The uncertainty 
on the integrated cross-section is approximately the same since 
this correction is a fully correlated yield adjustment. 

Table |2] summarizes these quantities for each of the var- 
ious error sources considered with a typical total uncertainty 
of 9.8% on the double-differential cross-section values being 
reported and a 4.9% uncertainty on the total integrated cross- 
section. 



7.2 Parametrization of pion production data 



Sanford and Wang [161 have developed an empirical parametriza- 
tion for describing the production cross-sections of mesons in 
proton-nucleus interactions. This parametrization has the func- 
tional form: 



d^o-(pH-A 



X) 



dpdf2 



(p,e)=exp[AK^(l- 



Phea 



(12) 
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Error Category 


Comment on method for estimating error 


(JJiff (%) 


S[nt (%) 


Statistical Errors: 








Be target statistics 


statistical error 


4.2 


0.6 


Empty target subtraction 


statistical error 


4.6 


0.6 


Sub-total 




6.3 


0.8 


Track yield corrections: 








Reconstruction efficiency 


Stat ot tracking efficiency computation sample 


1.3 


0.8 


Pion, proton absorption 


10% uncertainty on n, p absorption rates 


3.6 


3.7 


Tertiary subtraction 


50% uncertainty on tertiary production rate 


1.8 


1.8 


Empty target subtraction 


5% uncertainty on empty target subtraction normalization 


1.3 


1.2 


5>uD-totai 




4.0 




Particle Identification: 








Electron veto 


Stat of electron veto efficiency computation sample 


u.z 


<U. i 


Pion, proton ID correction 


analytical propagation of errors from parameterized 
PID detector response functions as described in 1111 


0.4 


0.1 


Sub-total 




0.5 


0.1 


Momentum reconstruction: 








Momentum scale 


2% uncertainty on absolute momentum scale 


3.6 


0.1 


Momentum resolution 


different hadronic generators used to generate correction 


3.4 


1.0 


Sub-total 




5.2 


1.0 


Overall normalization: 


targeting eff., fully correlated recon. and PID contributions 


2.0 


2.0 


Total 


9.8 


4.9 



Table 2. Summary of the uncertainties affecting the double-differential cross-section (5diff) and integrated cross-section (Sint) measurements 
(defined in text). 



where: 



Ci - C3- 



Ce9{p- CyPbeam COs'' 



(13) 



and X denotes any system of other particles in the final state, 
Pbeam '^^e proton beam momentum in GeV /c, p and 6 are 
the 7r+ momentum and angle in units of GeV /c and radians, re- 
spectively, cPa/ {dpdff) is expressed in units of mb/(GeV/c)/sr, 
dfi = 2tt (i(cos 9), and the parameters ci, . . . , cs are obtained 
from fits to meson production data. 

The parameter ci is an overall normalization factor, the four 
parameters C2, C3, C4, C5 describe the momentum distribution of 
the secondary pions in the forward direction, and the three pa- 
rameters cq,cj, cg describe the corrections to the pion momen- 
tum distribution for pion production angles that are different 
from zero. 

The 7r+ production data reported here have been fitted to 
the empirical Sanford-Wang formula. In the minimization 
procedure, seven out of these eight parameters were allowed to 
vary. The parameter C5 was fixed to the conventional value C5 = 
C4, since the cross-section dependence on the proton beam mo- 
mentum cannot be addressed by the present HARP data-set, 
which includes exclusively measurements taken at p^eam = 
8.9 GeV/c. In the minimization, the full error matrix was 
used. The goodness-of-fit of the Sanford-Wang parametriza- 
tion hypothesis for the HARP results can be assessed by con- 
sidering the best-fit value of x^jj^ = 248 for 7 1 degrees of 
freedom, indicating a very poor fit quality. In particular, inspec- 
tion of the HARP inclusive pion production double-differential 
cross-section, and resulting Sanford-Wang parametrization. 



points to a description of the ratio g{9) of the pion momen- 
tum distribution at 6' ^ with respect to the 9 = pion mo- 
mentum distribution that is more complicated than what can 
be accommodated within the Sanford-Wang formula, where 
this ratio is given by 17(6') = exp[—ce9{p — pc)], with pc = 
CTPbeam COS^« 9. 

Given the poor description of this HARP pion production 
data-set in terms of the original Sanford-Wang parametrization, 
we explored alternative functional forms. We found a signif- 
icantly better representation of the data by adopting a simple 
generalization of the Sanford-Wang formula, obtained by in- 
troducing one extra-parameter cg 7^ for the description of 
the angular dependence of the pion momentum distribution, ac- 



cording to .g(6l) = (1+p/Pbeam) 



C9S(p-Pc) 



exp[-C69{p~Pc)]- 



Overall, we use the following parametrization for the inclusive 
7r+ production double-differential cross-section: 



dpdil Pbeam 



x(l 



p 



Pbeam 



(14) 



where the argument A in the exponent is given by Eq. [T3] We 

obtain in this case a best-fit value of xt,- = 117 for 70 

-^^min 

degrees of freedom. 

Concerning the parameters estimation, the best-fit values of 
the extended Sanford-Wang parameter set discussed above are 
reported in Table[3] together with their eiTors. The fit parameter 
eiTors are estimated by requiring Ay^ = x^ ~ X^in ^ 9.30, 
corresponding to the 68.27% confidence level region for eight 
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Paranictcr 




Cl 


(5.13 ±0.41) 


C2 


(1.87 ±0.52) 


C3 


(6.67 ± 1.69) 


C4 — C5 


(1.56 ±0.55) 


C6 


(1.19 ±0.18) ■ 10^ 


c? 


(1.73 ±0.31) ■ 10"^ 


C8 


(1.98 ±0.69) ■ 10^ 


cg 


(1.60 ±0.44) ■ 10^ 



Table 3. Extended Sanford-Wang parameters and errors obtained by fitting the dataset. The errors refer to the 68.27% confidence level for 
eight parameters (Ax^ = 9.30). 



Parameter 


Cl 


C2 


C3 


C4 = C5 


C6 


C7 


C8 




Cl 


1.000 
















C2 


0.341 


1.000 














C3 


0.099 


0.562 


1.000 












C4 = C5 


-0.696 


-0.730 


0.004 


1.000 










C6 


-0.309 


0.288 


0.735 


0.358 


1.000 








C7 


-0.609 


0.066 


-0.221 


0.030 


0.005 


1.000 






C8 


-0.170 


-0.030 


-0.270 


-0.173 


-0.433 


0.672 


1.000 




C9 


-0.250 


0.270 


0.819 


0.368 


0.973 


-0.060 


-0.393 


1.000 



Table 4. Correlation coefficients among the extended Sanford-Wang parameters, obtained by fitting the data. 



variable parameters. Significant correlations among fit param- 
eters are found, as shown by the correlation matrix given in 
Table m 

The HARP cross-section measurement is compared to the 
best-fit parametrization of Eqs.[T4]and[T3] and Table|3] in Figs. [17] 
and [18] Also by looking at Fig.[T7] one can qualitatively con- 
clude that the proposed generalization of the Sanford-Wang 
parametrization provides a reasonable description of the data 
spectral features over the entire pion phase space measured. 
On the other hand, as already noted in 1 10], we remark that the 
goodness-of-fit depends on the correlations among the HARP 
cross-section uncertainties in different {p, 9) bins, and there- 
fore cannot be inferred solely from Fig. [17] 

We defer to a later publication, which should include a 
more comprehensive study of 7r+ production at various beam 
momenta and from various nuclear targets, a more complete 
discussion on the adequacy of parametrization-driven models 
such as Sanford-Wang (or simple modifications of) to describe 
HARP hadron production data. 



8 Relevance of HARP beryllium results for 
neutrino experiments 

The Booster neutrino beam at the Fermi National Acceler- 
ator Laboratory in Batavia, Illinois is created from the decay 
of charged mesons passing through a 50 m open decay region. 
These mesons are produced when 8.9 GeV/c momentum pro- 
tons are impinged upon a 71 cm long (1.7A) by 1 cm diame- 
ter beryllium target located at the upstream end of a magnetic 
focusing horn. This neutrino beam has been used by the Mini- 
BooNE experiment since September, 2002 and will be used by 
the SciBooNE experiment starting in summer, 2007. 



The MiniBooNE (E898) experiment at Fermilab ^ was 
designed to address the yet unconfirmed oscillation signal re- 
ported by the LSND collaboration |6|. MiniBooNE has been 
searching for the appearance of electron neutrinos in a beam 
that is predominantly muon flavor with an L/E similar to LSND 
but with substantially differing systematics. Additionally, the 
MiniBooNE detector can be used to make neutrino interac- 
tion cross-section measurements for both charged-current and 
neutral-current processes. An important systematic for these 
analyses arises from the prediction of the fluxes of different 
neutrino flavors at the MiniBooNE detector. For the ly^ appear- 
ance search the effect of the normalization uncertainty on 7r+ 
production is largely reduced by a constraint provided by the 
in situ measurement of charged-current quasi-elastic events. 
Neutrino cross-section measurements at MiniBooNE, however, 
will directly benefit from reductions in neutrino flux normaliza- 
tion uncertainties enabled by these hadron data. 

The neutrino flux prediction at the MiniBooNE detector 
is generated using a Monte Carlo simulation implemented in 
Geant4 IJSJ . Primary meson production rates are presently de- 
termined by fitting the empirical parametrization of Sanford 
and Wang [16] to production data in the relevant region. The 
results presented here, being for protons at exactly the Booster 
beam energy, are a critical addition to the global Sanford-Wang 
parametrization fits(3 

A flux prediction based on these data has been used by the 
MiniBooNE collaboration in their search for ly^ appearance at 
Am^ ^ 1 eV^ |5J. The systematic uncertainty on their pre- 



Cross-section data from the Brookhaven E910 experiment are also 
used to provide an additional constraint to the Sanford-Wang formula 
at angles larger than 210 mrad. These data are from 6.4 and 12.3 
GeV/c proton beams on a beryllium target. 
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Fig. 19. Fractional uncertainty (in percent) on the double-differential 
pion production cross-section measured, as a function of pion momen- 
tum and angle. The total uncertainty is shown by the thick black his- 
tograms, and individual contributions from the error categories given 
in Tab. |2] are also shown. Statistical, track yield corrections, and mo- 
mentum reconstruction uncertainties are shown as thin solid, dashed, 
and dotted histograms, respectively; the overall normalization uncer- 
tainty is not shown, and the particle identification uncertainty contri- 
bution lies below 2% for all pion momenta and angles. 



diction of Vr from 7r+ 



as a background to the 




Fig. 20. Predicted muon neutrino flux at the MiniBooNE detector 
from a Geant4 simulation of the Booster Neutrino Beam at Fermi- 
lab based on the parametrization of the HARP n'^ production cross- 
section (Eqs. [12] and [13] and Table [3j. The solid curve is the total 
muon neutrino flux, while the dashed curve is the part of the flux 
coming from the decay of created in proton-beryllium collisions. 
The primary production of positive pions is based on a parametriza- 
tion of the HARP tt'^ cross-section measurements presented in this 
paper and represents the predominant source of at MiniBooNE. 
The dotted histogram shows the part of the flux coming from 
the decay of 7r"'"s that are within the kinematic boundaries of the 
measurement presented here, 0.75 GeV/c < p-n < 6.5 GeV/c and 
0.030 rad < 6^ < 0.210 rad. 
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oscillation analysis is 8% [5J with only a 3% contribution from 
the 7r+ production model based on these data. In contrast, the 
uncertainty on the prediction of Ve from kaon decays is 35% 
|5], the largest contribution (25%) coming from the production 
model of kaons in pH-Be interactions, also based on parameter- 
izations of available cross-section data. 

Using the complete MiniBooNE beam Monte Carlo we can 
illustrate the direct impact of the HARP data on the MiniBooNE 
flux predictions. The dominant channel leading to a muon neu- 
trino in the detector is p + Be ^ 7r+ — ^ v^. Fig.l20lshows the 
total z/p flux (solid) according to the simulation as well as the 
part coming directly from the sequence listed above (dashed). 
The curves are obtained from a Geant4 simulation of the Booster 
Neutrino Beam at Fermilab based on the parametrization of the 
HARP 7r+ production cross-section given in Eqs. [12] and [13] 
and Table[3] Figure|2T]shows the kinematic distribution of tt+'s 
which result in a t'p in the MiniBooNE detector The box out- 
lines the kinematic range of the measurements described in this 
paper. The simulation indicates that >80% of the relevant pions 
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Fig. 21. p-9 distribution of 7r"'"'s which decay to a muon neutrino 
passing through the MiniBooNE detector according to a Monte Carlo 
simulation. The box marks the kinematic region of the measurement 
presented here and contains 81.1% of the pions contributing to the 
neutrino flux at MiniBooNE. 



come from within this region. The neutrinos produced by the 
subset of pion phase-space directly covered by this measure- 
ment are shown by the dotted histogram in Fig. [20] While the 
coverage of these data is being displayed for the MiniBooNE 
detector, a similar coverage is expected for the SciBooNE de- 
tector in located in the same neutrino beam at Fermilab. 
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9 Summary and conclusions 

In this paper we have presented a measurement of the double- 
differential production cross-section of positive pions in the 
collision of 8.9 GeV/c protons with a beryllium target. The 
data have been reported in 78 bins of pion momentum and 
angle in the kinematic range from 0.75 GeV/c < p„ < 6.5 
GeV/c and 0.030 rad < 0^ < 0.210 rad. A systematic error 
analysis has been performed yielding an average point-to-point 
error of 9.8% (statistical + systematic) and an uncertainty on 
the total integrated cross-section of 4.9%. Further, the data have 
been fitted to a modified form of the empirical parameterization 
of Sanford and Wang and the resulting parameters provided. 

These production data have direct relevance for the predic- 
tion of a flux for MiniBooNE, an experiment searching for 
i^fi t^e oscillations using the Booster neutrino beam line at 
Fermi National Accelerator Laboratory, and SciBooNE, an ex- 
periment designed to measure cross-sections in the 1 GeV 
neutrino energy region using the same beam. Final flux pre- 
dictions for these experiments will be based on the results pre- 
sented here and published elsewhere by the MiniBooNE and 
SciBooNE coUaborations. 
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Table 5. HARP results for the double-differential tt"*" production 
cross-section in the laboratory system, d^a^ /{dpdO). Each row 
refers to a different (Pjjjjn < p < pmax, ^'jjjj^ < < 9max) bin, 
where p and 6 are the pion momentum and polar angle, respectively. 
The central value as well as the square-root of the diagonal elements 
of the covariance matrix are given. 



